Efficient algorithms for rigid body integration using optimized splitting methods and exact free rotational motion.

In this note we present molecular dynamics integration schemes that combine optimized splitting and gradient methods with exact free rotational motion for rigid body systems and discuss their relative merits. The algorithms analyzed here are based on symplectic, timereversible schemes that conserve all relevant constants of the motion. It is demonstrated that although the algorithms differ in their stability due to truncation errors associated with limited numerical precision, the optimized splitting methods can outperform the commonly-used velocity Verlet scheme at a level of precision typical of most simulations in which dynamical quantities are of interest. Useful guidelines for choosing the best integration scheme for a given level of accuracy and stability are provided. Hamiltonian splitting methods are an established technique to derive stable and accurate integration schemes in molecular dynamics. The strategy of these methods is to split the Hamiltonian of the system into parts whose evolution can be solved exactly. Using the Campbell-BakerHausdorff formula, splitting algorithms can be presented as products of exactly solvable propagation steps, involving more factors for higher-order schemes. The resulting algorithms can be optimized by adjusting the form of the splitting to minimize error estimates. Recently, secondand fourth-order symplectic integration schemes for simulations of rigid body motion, based on the exact solution for the full kinetic (free) propagator, have been proposed. While this exact solution involves elliptic functions, elliptic integrals and theta functions, there exist efficient numerical routines to compute elliptic functions, and the computation of elliptic integrals and theta functions can be implemented efficiently or avoided altogether using a recursive method. Employing the exact free rotational motion, the resulting splitting method leads to demonstrably more accurate dynamics for systems in which free motion is important. Furthermore, using the exact kinetic propagator, any splitting scheme for integrating the dynamics of point particles can be transferred to rigid systems. Here we analyze the combination of the exact kinetic propagator and optimized splitting and gradient-like approaches. For a system of rigid bodies, a phase space point Γ is specified by a center of mass position qi, an attitude matrix Si, and translational and angular momenta pi and li for each particle i of mass mi. Given the Hamiltonian H = T +V , where T and V are the kinetic and potential energies, respectively, the time evolution of the point Γ in phase space is governed by Γ̇ = {H,Γ} = {T,Γ}+{V,Γ}, in which {, } denotes the Poisson bracket. Henceforth, the operators {T, .} and {V, .} will be designated asA and B, respectively. Defining L = A+B, the solution of the equations of motion is formally given by Γ(t) = eΓ(0). While the various possible splitting schemes can be assigned a theoretical efficiency, the relative efficiency of real simulations can be somewhat different. Nonetheless, the estimates are useful to eliminate the least efficient variants. Based on our studies of second and fourth order methods, the most efficient integration schemes can be formulated using the following generic form of the splitting algorithm for a single time step of size h: